Theme Song
# install.packages("remotes")
require('BBmisc')
#remotes::install_github("rstudio/sass")
lib('sass')
## sass
## TRUE
/* https://stackoverflow.com/a/66029010/3806250 */
h1 { color: #002C54; }
h2 { color: #2F496E; }
h3 { color: #375E97; }
/* ----------------------------------------------------------------- */
/* https://gist.github.com/himynameisdave/c7a7ed14500d29e58149#file-broken-gradient-animation-less */
.hover01 {
/* color: #FFD64D; */
background: linear-gradient(155deg, #EDAE01 0%, #FFEB94 100%);
transition: all 0.45s;
&:hover{
background: linear-gradient(155deg, #EDAE01 20%, #FFEB94 80%);
}
}
.hover02 {
color: #FFD64D;
background: linear-gradient(155deg, #002C54 0%, #4CB5F5 100%);
transition: all 0.45s;
&:hover{
background: linear-gradient(155deg, #002C54 20%, #4CB5F5 80%);
}
}
.hover03 {
color: #FFD64D;
background: linear-gradient(155deg, #A10115 0%, #FF3C5C 100%);
transition: all 0.45s;
&:hover{
background: linear-gradient(155deg, #A10115 20%, #FF3C5C 80%);
}
}
## https://stackoverflow.com/a/36846793/3806250
options(width = 999)
knitr::opts_chunk$set(class.source = 'hover01', class.output = 'hover02', class.error = 'hover03')
課題をすぐに提出してください
課題の提出期限は、5月31日 15:59 JSTですが、可能であれば1日か2日早く提出してください。 早い段階で提出すると、他の受講生のレビューを時間内に得る可能性が高くなります。
The R dataset faithful contains data on waiting time between eruptions (the column named waiting) and the duration of the eruption (the column named eruptions) for the famous Old Faithful geyser in Yellowstone National Park, Wyoming, USA.
In this case, you are asked to modify the MCMC algorithm provided in “Sample code for density estimation problems” (as opposed to the EM algorithm you used in the previous peer assignment) to provide a (Bayesian) density estimate the marginal distribution of the duration of the eruptions using a location-and-scale mixture of 2 univariate Gaussian distributions (as opposed to the location mixture of 6 univariate Gaussian distributions that we used for the galaxies dataset). Assume that the priors are \(\omega∼Beta(1,1)\), \(\mu_k∼Normal(η,\tau^2)\) and \(1/σ^2_k∼Gamma(d,q)\), where \(\eta\), \(\tau^2\), \(d\) and \(q\) are selected using an empirical Bayes approach similar to the one we used in “Sample code for density estimation problems”.
Reviewers will check whether the code has been modified correctly, and whether the density estimate you generate appears correct. Please remember that you are being asked to use a location-and-scale mixture to generate the density estimate, so the “Sample code for density estimation problems” cannot be used directly and requires some modification. Before submitting your answer, it might be useful to compare the density estimate generated by your algorithm against a kernel density estimate generated by the R function density(), and agains the answer to the previous peer assignment.
if(!suppressPackageStartupMessages(require('BBmisc'))) {
install.packages('BBmisc', dependencies = TRUE, INSTALL_opts = '--no-lock')
}
suppressPackageStartupMessages(require('BBmisc'))
# suppressPackageStartupMessages(require('rmsfuns'))
pkgs <- c('devtools', 'knitr', 'kableExtra', 'tidyr',
'readr', 'lubridate', 'data.table', 'reprex',
'timetk', 'plyr', 'dplyr', 'stringr', 'magrittr',
'tdplyr', 'tidyverse', 'formattable',
'echarts4r', 'paletteer')
suppressAll(lib(pkgs))
# load_pkg(pkgs)
## Set the timezone but not change the datetime
Sys.setenv(TZ = 'Asia/Tokyo')
## options(knitr.table.format = 'html') will set all kableExtra tables to be 'html', otherwise need to set the parameter on every single table.
options(warn = -1, knitr.table.format = 'html')#, digits.secs = 6)
## https://stackoverflow.com/questions/39417003/long-vectors-not-supported-yet-abnor-in-rmd-but-not-in-r-script
knitr::opts_chunk$set(message = FALSE, warning = FALSE)#,
#cache = TRUE, cache.lazy = FALSE)
rm(pkgs)
dat <- fread('data/fuses.csv') %>% as.matrix
head(dat)
## V1
## [1,] 1.062163
## [2,] 1.284631
## [3,] 2.733352
## [4,] 2.284062
## [5,] 3.289444
## [6,] 3.113181
Source : 400 x 1
Is the sampler for \(c\) correct?
The full conditional for \(c_i\) is given by,
\[\begin{align} \Pr(c_i = 1|...) = \frac{\omega\frac{1}{\sqrt[]{2\pi\sigma_1}}exp\left\{-\frac{1}{2\sigma^2_1}(x_i-\mu_1)^2\right\}}{\omega\frac{1}{\sqrt[]{2\pi\sigma_1}}exp\left\{-\frac{1}{2\sigma^2_1}(x_i-\mu_1)^2\right\}+(1-\omega)\frac{1}{\sqrt[]{2\pi\sigma_2}}exp\left\{-\frac{1}{2\sigma^2_2}(x_i-\mu_2)^2\right\}} \end{align}\]
(This assumes that the mixture weights are parameterized as
for(i in 1:n){
v = rep(0,2)
v[1] = log(w) + dnorm(x[i], mu[1], sigma[1], log=TRUE) #Compute the log of the weights
v[2] = log(1-w) + dnorm(x[i], mu[2], sigma[2], log=TRUE) #Compute the log of the weights
v = exp(v - max(v))/sum(exp(v - max(v)))
cc[i] = sample(1:2, 1, replace=TRUE, prob=v)
}
## Error in 1:n: NA/NaN argument
Is the sampler for \(\omega\) correct?
The full conditional is simply
\[ \omega \mid \cdots \mbox{Beta} \left(\alpha_1 + \sum_{i=1}^{n} \mathbf{1}(c_i=1), \alpha_2 + \sum_{i=1}^{n} \mathbf{1}(c_i=2)\right) \]
The following one line of code implements it:
w = rbeta(1, aa[1] + sum(cc==1), aa[2] + sum(cc==2))
## Error in rbeta(1, aa[1] + sum(cc == 1), aa[2] + sum(cc == 2)): object 'aa' not found
Are the sampler for \(\mu\) correct?
\[ \mu_k \mid \cdots \mbox{Normal}\left( \left[ \frac{n_k}{\sigma_k^2} + \frac{1}{\tau^2} \right]^{-1}\left[ \frac{\sum_{i:c_i=k} x_i}{\sigma_k^2} + \frac{\eta}{\tau^2} \right], \left[ \frac{n_k}{\sigma_k^2} + \frac{1}{\tau^2} \right]^{-1}\right) \]
The following code implements the sampler:
# Sample the means
for(k in 1:2){
nk = sum(cc==k)
xsumk = sum(x[cc==k])
tau2.hat = 1/(nk/sigma[k]^2 + 1/tau^2)
mu.hat = tau2.hat*(xsumk/sigma[k]^2 + eta/tau^2)
mu[k] = rnorm(1, mu.hat, sqrt(tau2.hat))
}
## Error in eval(expr, envir, enclos): object 'cc' not found
Remember, however, that there are various ways in which the implementation could proceed.
Are the samplers for \(\sigma_k\) correct?
The full conditional takes the form:
\[ 1/\sigma^2_k \mid \cdots \mbox{Gamma}\left(d + \frac{n_k}{2} , q + \frac{1}{2} \sum_{i:c_k=k} (x_i - \mu_k)^2 \right) \]
The following code implements the sampler:
# Sample the variances
for(k in 1:2){
nk = sum(cc==k)
dd.star = dd + nk/2
qq.star = qq + sum((x[cc==k] - mu[k])^2)/2
sigma[k] = sqrt(1/rgamma(1, dd.star, qq.star))
}
## Error in eval(expr, envir, enclos): object 'cc' not found
Remember, however, that there are various ways in which the implementation could proceed.
Does the code to generate the density estimate appear correct?
This is an example that, as before, assumes that the weights are parameterized as \(\omega\) and \(1- \omega\).
## Plot Bayesian estimate with pointwise credible bands
density.mcmc.lq = apply(density.mcmc, 2, quantile, 0.025)
## Error in apply(density.mcmc, 2, quantile, 0.025): object 'density.mcmc' not found
density.mcmc.uq = apply(density.mcmc, 2, quantile, 0.975)
## Error in apply(density.mcmc, 2, quantile, 0.975): object 'density.mcmc' not found
plot(xx, density.mcmc.m, type="n",ylim=c(0,max(density.mcmc.uq)),xlab="Eruptions", ylab="Density")
## Error in plot(xx, density.mcmc.m, type = "n", ylim = c(0, max(density.mcmc.uq)), : object 'xx' not found
polygon(c(xx,rev(xx)), c(density.mcmc.lq, rev(density.mcmc.uq)), col="grey", border="grey")
## Error in xy.coords(x, y, setLab = FALSE): object 'density.mcmc.lq' not found
lines(xx, density.mcmc.m, col="black", lwd=2)
## Error in lines(xx, density.mcmc.m, col = "black", lwd = 2): object 'xx' not found
points(x, rep(0,n))
## Error in points(x, rep(0, n)): object 'x' not found
Note that this will generate not only a point estimate, but also pointwise credible intervals.
Does the graph appear to be correct?
The density estimate should look similar to the Figure below. Note that the variance of both components is clearly different.
The full code for the MCMC algorithm follows:
x = faithful$eruptions
n = length(x)
plot(density(x))
points(x,rep(0,n))
## Priors set up using an "empirical Bayes" approach
aa = rep(1,2)
eta = mean(x)
tau = sqrt(var(x))
dd = 2
qq = var(x)/2
## Initialize the parameters
w = 1/2
mu = rnorm(2, mean(x), sd(x))
sigma = rep(sd(x)/2,2)
cc = sample(1:2, n, replace=T, prob=c(1/2,1/2))
## Number of iterations of the sampler
rrr = 12000
burn = 2000
## Storing the samples
cc.out = array(0, dim=c(rrr, n))
w.out = rep(0, rrr)
mu.out = array(0, dim=c(rrr, 2))
sigma.out = array(0, dim=c(rrr, 2))
logpost = rep(0, rrr)
for(s in 1:rrr){
# Sample the indicators
for(i in 1:n){
v = rep(0,2)
v[1] = log(w) + dnorm(x[i], mu[1], sigma[1], log=TRUE) #Compute the log of the weights
v[2] = log(1-w) + dnorm(x[i], mu[2], sigma[2], log=TRUE) #Compute the log of the weights
v = exp(v - max(v))/sum(exp(v - max(v)))
cc[i] = sample(1:2, 1, replace=TRUE, prob=v)
}
# Sample the weights
w = rbeta(1, aa[1] + sum(cc==1), aa[2] + sum(cc==2))
# Sample the parameters of the components
for(k in 1:2){
# Sample the means
nk = sum(cc==k)
xsumk = sum(x[cc==k])
tau2.hat = 1/(nk/sigma[k]^2 + 1/tau^2)
mu.hat = tau2.hat*(xsumk/sigma[k]^2 + eta/tau^2)
mu[k] = rnorm(1, mu.hat, sqrt(tau2.hat))
# Sample the variances
dd.star = dd + nk/2
qq.star = qq + sum((x[cc==k] - mu[k])^2)/2
sigma[k] = sqrt(1/rgamma(1, dd.star, qq.star))
}
# Store samples
cc.out[s,] = cc
w.out[s] = w
mu.out[s,] = mu
sigma.out[s,] = sigma
logpost[s] = 0
for(i in 1:n){
if(cc[i] == 1){
logpost[s] = logpost[s] + log(w) + dnorm(x[i], mu[1], sigma[1], log=TRUE)
}else{
logpost[s] = logpost[s] + log(1-w) + dnorm(x[i], mu[2], sigma[2], log=TRUE)
}
}
logpost[s] = logpost[s] + dbeta(w, aa[1], aa[2], log=TRUE)
for(k in 1:2){
logpost[s] = logpost[s] + dnorm(mu[k], eta, tau, log = T) + dgamma(1/sigma[k]^2, dd, qq)/sigma[k]^4
}
if(s/500==floor(s/500)){
print(paste("s =",s))
}
}
## [1] "s = 500"
## [1] "s = 1000"
## [1] "s = 1500"
## [1] "s = 2000"
## [1] "s = 2500"
## [1] "s = 3000"
## [1] "s = 3500"
## [1] "s = 4000"
## [1] "s = 4500"
## [1] "s = 5000"
## [1] "s = 5500"
## [1] "s = 6000"
## [1] "s = 6500"
## [1] "s = 7000"
## [1] "s = 7500"
## [1] "s = 8000"
## [1] "s = 8500"
## [1] "s = 9000"
## [1] "s = 9500"
## [1] "s = 10000"
## [1] "s = 10500"
## [1] "s = 11000"
## [1] "s = 11500"
## [1] "s = 12000"
## Compute the samples of the density over a dense grid
xx = seq(0,7,length=150)
density.mcmc = array(0, dim=c(rrr-burn,length(xx)))
for(s in 1:(rrr-burn)){
density.mcmc[s,] = density.mcmc[s,] + w.out[s+burn]*dnorm(xx,mu.out[s+burn,1],sigma.out[s+burn,1]) +
(1-w.out[s+burn])*dnorm(xx,mu.out[s+burn,2],sigma.out[s+burn,2])
}
density.mcmc.m = apply(density.mcmc , 2, mean)
## Plot Bayesian estimate with pointwise credible bands
density.mcmc.lq = apply(density.mcmc, 2, quantile, 0.025)
density.mcmc.uq = apply(density.mcmc, 2, quantile, 0.975)
plot(xx, density.mcmc.m, type="n",ylim=c(0,max(density.mcmc.uq)),xlab="Eruptions", ylab="Density")
polygon(c(xx,rev(xx)), c(density.mcmc.lq, rev(density.mcmc.uq)), col="grey", border="grey")
lines(xx, density.mcmc.m, col="black", lwd=2)
points(x, rep(0,n))
prompt caption-text color-secondary-text
### Get a "Bayesian" kernel density estimator based on the same location mixture of 6 normals
## Priors set up using an "empirical Bayes" approach
aa = rep(1,KK)
## Error in eval(expr, envir, enclos): object 'KK' not found
eta = mean(x)
tau = sqrt(var(x))
dd = 2
qq = var(x)/KK
## Error in eval(expr, envir, enclos): object 'KK' not found
mu_0 = rnorm(KK, eta, tau)
## Error in rnorm(KK, eta, tau): object 'KK' not found
sigma_0 = rinvgamma(KK, dd, qq)
## Error in rinvgamma(KK, dd, qq): could not find function "rinvgamma"
## Initialize the parameters
w = rep(1,KK)/KK
## Error in eval(expr, envir, enclos): object 'KK' not found
mu = rnorm(KK, mean(x), sd(x))
## Error in rnorm(KK, mean(x), sd(x)): object 'KK' not found
sigma = sd(x)/KK
## Error in eval(expr, envir, enclos): object 'KK' not found
cc = sample(1:KK, n, replace=T, prob=w)
## Error in sample(1:KK, n, replace = T, prob = w): object 'KK' not found
## Number of iterations of the sampler
rrr = 12000
burn = 2000
## Storing the samples
cc.out = array(0, dim=c(rrr, n))
w.out = array(0, dim=c(rrr, KK))
## Error in array(0, dim = c(rrr, KK)): object 'KK' not found
mu.out = array(0, dim=c(rrr, KK))
## Error in array(0, dim = c(rrr, KK)): object 'KK' not found
sigma.out = rep(0, rrr)
logpost = rep(0, rrr)
for(s in 1:rrr){
# Sample the indicators
for(i in 1:n){
v = rep(0, KK)
for(k in 1:KK){
v[k] = log(w[k]) + dnorm(x[i], mu[k], sigma, log = TRUE) #Compute the log of the weights
}
v = exp(v - max(v))/sum(exp(v - max(v)))
cc[i] = sample(1:KK, 1, replace = TRUE, prob = v)
}
# Sample the weights
w = as.vector(rdirichlet(1, aa + tabulate(cc, nbins=KK)))
# Sample the means
for(k in 1:KK){
nk = sum(cc==k)
xsumk = sum(x[cc==k])
tau2.hat = 1/(nk/sigma^2 + 1/sigma_0[k]^2)
mu.hat = tau2.hat*(xsumk/sigma^2 + mu_0[k]/sigma_0[k]^2)
mu[k] = rnorm(1, mu.hat, sqrt(tau2.hat))
}
# Sample the variances
dd.star = dd + n/2
qq.star = qq + sum((x - mu[cc])^2)/2
sigma = sqrt(1/rgamma(1, dd.star, qq.star))
# Store samples
cc.out[s,] = cc
w.out[s,] = w
mu.out[s,] = mu
sigma.out[s] = sigma
for(i in 1:n){
logpost[s] = logpost[s] + log(w[cc[i]]) + dnorm(x[i], mu[cc[i]], sigma, log = TRUE)
}
logpost[s] = logpost[s] + log(ddirichlet(w, aa))
for(k in 1:KK){
logpost[s] = logpost[s] + dnorm(mu[k], eta, tau, log = TRUE)
}
logpost[s] = logpost[s] + dgamma(1/sigma^2, dd, qq, log = TRUE) - 4*log(sigma)
if(s/500==floor(s/500)){
print(paste("s =",s))
}
}
## Error in eval(expr, envir, enclos): object 'KK' not found
Provide code to generate the density estimate on a grid consisting of 150 points in the interval [0,7].
xx = seq(0,7,length=150)
nxx = length(xx)
## Compute the samples of the density over a dense grid
density.mcmc = array(0, dim = c(rrr-burn, length(xx)))
for(s in 1:(rrr-burn)){
for(k in 1:KK){
density.mcmc[s,] = density.mcmc[s,] + w.out[s+burn,k]*dnorm(xx,mu.out[s+burn,k],sigma.out[s+burn])
}
}
## Error in eval(expr, envir, enclos): object 'KK' not found
density.mcmc.m = apply(density.mcmc , 2, mean)
Provide the a graph of the density estimate on the interval [0,7].
## Plot Bayesian estimate with pointwise credible bands along with kernel density estimate and frequentist point estimate
colscale = c("black", "blue", "red")
yy = density(x)
density.mcmc.lq = apply(density.mcmc, 2, quantile, 0.025)
density.mcmc.uq = apply(density.mcmc, 2, quantile, 0.975)
plot(xx, density.mcmc.m, type="n",ylim=c(0,max(density.mcmc.uq)),xlab="Eruption Duration", ylab="Density")
polygon(c(xx,rev(xx)), c(density.mcmc.lq, rev(density.mcmc.uq)), col="grey", border="grey")
lines(xx, density.mcmc.m, col=colscale[1], lwd=2)
#lines(xx, density.EM, col=colscale[2], lty=2, lwd=2)
lines(yy, col=colscale[3], lty=3, lwd=2)
points(x, rep(0,n))
legend(5, 0.45, c("KDE","EM","MCMC"), col=colscale[c(3,2,1)], lty=c(3,2,1), lwd=2, bty="n")
Is the sampler for \(c\) correct?
The full conditional for \(c_i\) is given by,
\[ \Pr(c_i = 1|...) = \frac{\omega\frac{1}{\sqrt[]{2\pi\sigma_1}}exp\left\{-\frac{1}{2\sigma^2_1}(x_i-\mu_1)^2\right\}}{\omega\frac{1}{\sqrt[]{2\pi\sigma_1}}exp\left\{-\frac{1}{2\sigma^2_1}(x_i-\mu_1)^2\right\}+(1-\omega)\frac{1}{\sqrt[]{2\pi\sigma_2}}exp\left\{-\frac{1}{2\sigma^2_2}(x_i-\mu_2)^2\right\}} \]
(This assumes that the mixture weights are parameterized as
for(i in 1:n){
v = rep(0,2)
v[1] = log(w) + dnorm(x[i], mu[1], sigma[1], log=TRUE) #Compute the log of the weights
v[2] = log(1-w) + dnorm(x[i], mu[2], sigma[2], log=TRUE) #Compute the log of the weights
v = exp(v - max(v))/sum(exp(v - max(v)))
cc[i] = sample(1:2, 1, replace=TRUE, prob=v)
}
Is the sampler for \(\omega\) correct?
The full conditional is simply
\[ \omega \mid \cdots \mbox{Beta} \left(\alpha_1 + \sum_{i=1}^{n} \mathbf{1}(c_i=1), \alpha_2 + \sum_{i=1}^{n} \mathbf{1}(c_i=2)\right) \]
The following one line of code implements it:
w = rbeta(1, aa[1] + sum(cc==1), aa[2] + sum(cc==2))
Are the sampler for \(\mu\) correct?
\[ \mu_k \mid \cdots \mbox{Normal}\left( \left[ \frac{n_k}{\sigma_k^2} + \frac{1}{\tau^2} \right]^{-1}\left[ \frac{\sum_{i:c_i=k} x_i}{\sigma_k^2} + \frac{\eta}{\tau^2} \right], \left[ \frac{n_k}{\sigma_k^2} + \frac{1}{\tau^2} \right]^{-1}\right) \]
The following code implements the sampler:
# Sample the means
for(k in 1:2){
nk = sum(cc==k)
xsumk = sum(x[cc==k])
tau2.hat = 1/(nk/sigma[k]^2 + 1/tau^2)
mu.hat = tau2.hat*(xsumk/sigma[k]^2 + eta/tau^2)
mu[k] = rnorm(1, mu.hat, sqrt(tau2.hat))
}
Remember, however, that there are various ways in which the implementation could proceed.
Are the samplers for \(\sigma_k\) correct?
The full conditional takes the form:
\[ 1/\sigma^2_k \mid \cdots \mbox{Gamma}\left(d + \frac{n_k}{2} , q + \frac{1}{2} \sum_{i:c_k=k} (x_i - \mu_k)^2 \right) \]
The following code implements the sampler:
# Sample the variances
for(k in 1:2){
nk = sum(cc==k)
dd.star = dd + nk/2
qq.star = qq + sum((x[cc==k] - mu[k])^2)/2
sigma[k] = sqrt(1/rgamma(1, dd.star, qq.star))
}
Remember, however, that there are various ways in which the implementation could proceed.
Does the code to generate the density estimate appear correct?
This is an example that, as before, assumes that the weights are parameterized as \(\omega\) and \(1- \omega\).
## Plot Bayesian estimate with pointwise credible bands
density.mcmc.lq = apply(density.mcmc, 2, quantile, 0.025)
density.mcmc.uq = apply(density.mcmc, 2, quantile, 0.975)
plot(xx, density.mcmc.m, type="n",ylim=c(0,max(density.mcmc.uq)),xlab="Eruptions", ylab="Density")
polygon(c(xx,rev(xx)), c(density.mcmc.lq, rev(density.mcmc.uq)), col="grey", border="grey")
lines(xx, density.mcmc.m, col="black", lwd=2)
points(x, rep(0,n))
Note that this will generate not only a point estimate, but also pointwise credible intervals.
Does the graph appear to be correct?
The density estimate should look similar to the Figure below. Note that the variance of both components is clearly different.
The full code for the MCMC algorithm follows:
x = faithful$eruptions
n = length(x)
plot(density(x))
points(x,rep(0,n))
## Priors set up using an "empirical Bayes" approach
aa = rep(1,2)
eta = mean(x)
tau = sqrt(var(x))
dd = 2
qq = var(x)/2
## Initialize the parameters
w = 1/2
mu = rnorm(2, mean(x), sd(x))
sigma = rep(sd(x)/2,2)
cc = sample(1:2, n, replace=T, prob=c(1/2,1/2))
## Number of iterations of the sampler
rrr = 12000
burn = 2000
## Storing the samples
cc.out = array(0, dim=c(rrr, n))
w.out = rep(0, rrr)
mu.out = array(0, dim=c(rrr, 2))
sigma.out = array(0, dim=c(rrr, 2))
logpost = rep(0, rrr)
for(s in 1:rrr){
# Sample the indicators
for(i in 1:n){
v = rep(0,2)
v[1] = log(w) + dnorm(x[i], mu[1], sigma[1], log=TRUE) #Compute the log of the weights
v[2] = log(1-w) + dnorm(x[i], mu[2], sigma[2], log=TRUE) #Compute the log of the weights
v = exp(v - max(v))/sum(exp(v - max(v)))
cc[i] = sample(1:2, 1, replace=TRUE, prob=v)
}
# Sample the weights
w = rbeta(1, aa[1] + sum(cc==1), aa[2] + sum(cc==2))
# Sample the parameters of the components
for(k in 1:2){
# Sample the means
nk = sum(cc==k)
xsumk = sum(x[cc==k])
tau2.hat = 1/(nk/sigma[k]^2 + 1/tau^2)
mu.hat = tau2.hat*(xsumk/sigma[k]^2 + eta/tau^2)
mu[k] = rnorm(1, mu.hat, sqrt(tau2.hat))
# Sample the variances
dd.star = dd + nk/2
qq.star = qq + sum((x[cc==k] - mu[k])^2)/2
sigma[k] = sqrt(1/rgamma(1, dd.star, qq.star))
}
# Store samples
cc.out[s,] = cc
w.out[s] = w
mu.out[s,] = mu
sigma.out[s,] = sigma
logpost[s] = 0
for(i in 1:n){
if(cc[i] == 1){
logpost[s] = logpost[s] + log(w) + dnorm(x[i], mu[1], sigma[1], log=TRUE)
}else{
logpost[s] = logpost[s] + log(1-w) + dnorm(x[i], mu[2], sigma[2], log=TRUE)
}
}
logpost[s] = logpost[s] + dbeta(w, aa[1], aa[2], log=TRUE)
for(k in 1:2){
logpost[s] = logpost[s] + dnorm(mu[k], eta, tau, log = T) + dgamma(1/sigma[k]^2, dd, qq)/sigma[k]^4
}
if(s/500==floor(s/500)){
print(paste("s =",s))
}
}
## [1] "s = 500"
## [1] "s = 1000"
## [1] "s = 1500"
## [1] "s = 2000"
## [1] "s = 2500"
## [1] "s = 3000"
## [1] "s = 3500"
## [1] "s = 4000"
## [1] "s = 4500"
## [1] "s = 5000"
## [1] "s = 5500"
## [1] "s = 6000"
## [1] "s = 6500"
## [1] "s = 7000"
## [1] "s = 7500"
## [1] "s = 8000"
## [1] "s = 8500"
## [1] "s = 9000"
## [1] "s = 9500"
## [1] "s = 10000"
## [1] "s = 10500"
## [1] "s = 11000"
## [1] "s = 11500"
## [1] "s = 12000"
## Compute the samples of the density over a dense grid
xx = seq(0,7,length=150)
density.mcmc = array(0, dim=c(rrr-burn,length(xx)))
for(s in 1:(rrr-burn)){
density.mcmc[s,] = density.mcmc[s,] + w.out[s+burn]*dnorm(xx,mu.out[s+burn,1],sigma.out[s+burn,1]) +
(1-w.out[s+burn])*dnorm(xx,mu.out[s+burn,2],sigma.out[s+burn,2])
}
density.mcmc.m = apply(density.mcmc , 2, mean)
## Plot Bayesian estimate with pointwise credible bands
density.mcmc.lq = apply(density.mcmc, 2, quantile, 0.025)
density.mcmc.uq = apply(density.mcmc, 2, quantile, 0.975)
plot(xx, density.mcmc.m, type="n",ylim=c(0,max(density.mcmc.uq)),xlab="Eruptions", ylab="Density")
polygon(c(xx,rev(xx)), c(density.mcmc.lq, rev(density.mcmc.uq)), col="grey", border="grey")
lines(xx, density.mcmc.m, col="black", lwd=2)
points(x, rep(0,n))
Provide an MCMC algorithm to fit the two-component, location-and-scale mixture of Gaussians.
### Loading data and setting up global variables
library(MASS)
library(MCMCpack)
data(faithful)
KK = 2 # As asked
x = faithful[,1]
n = length(x)
set.seed(781209)
### Get a "Bayesian" kernel density estimator based on the same location mixture of 2 normals
## Priors set up using an "empirical Bayes" approach
aa = rep(1,KK)
eta = mean(x)
tau = sqrt(var(x))
dd = 2
qq = var(x)/KK
w = rep(1,KK)/KK
mu = rnorm(KK, mean(x), sd(x))
sigma = rep(sd(x)/KK, KK)
cc = sample(1:KK, n, replace=T, prob=w)
## Number of iterations of the sampler
rrr = 12000
burn = 2000
## Storing the samples
cc.out = array(0, dim=c(rrr, n))
w.out = array(0, dim=c(rrr, KK))
mu.out = array(0, dim=c(rrr, KK))
sigma.out = array(0, dim=c(rrr, KK))
logpost = rep(0, rrr)
for(s in 1:rrr){
# Sample the indicators
for(i in 1:n){
v = rep(0,KK)
for(k in 1:KK){
v[k] = log(w[k]) + dnorm(x[i], mu[k], sigma[k], log=TRUE) #Compute the log of the weights
}
v = exp(v - max(v))/sum(exp(v - max(v)))
cc[i] = sample(1:KK, 1, replace=TRUE, prob=v)
}
# Sample the weights
w = as.vector(rdirichlet(1, aa + tabulate(cc, nbins=KK)))
# Sample the means
for(k in 1:KK){
nk = sum(cc==k)
xsumk = sum(x[cc==k])
tau2.hat = 1/(nk/sigma[k]^2 + 1/tau^2)
mu.hat = tau2.hat*(xsumk/sigma[k]^2 + eta/tau^2)
mu[k] = rnorm(1, mu.hat, sqrt(tau2.hat))
}
# Sample the variances
for(k in 1:KK) {
dd.star = dd + sum(cc==k)/2
sumsq = 0
for(i in 1:n) {
if (cc[i]==k) {
sumsq = sumsq + ((x[i]-mu[k])^2)
}
}
qq.star = qq + sumsq/2
sigma[k] = sqrt(rinvgamma(1, shape = dd.star, scale = qq.star))
}
# Store samples
cc.out[s,] = cc
w.out[s,] = w
mu.out[s,] = mu
sigma.out[s,] = sigma
for(i in 1:n){
logpost[s] = logpost[s] + log(w[cc[i]]) + dnorm(x[i], mu[cc[i]], sigma[cc[i]], log=TRUE)
}
logpost[s] = logpost[s] + log(ddirichlet(w, aa))
for(k in 1:KK){
logpost[s] = logpost[s] + dnorm(mu[k], eta, tau, log=TRUE)
logpost[s] = logpost[s] + log(dinvgamma(sigma[k]^2, shape = dd, scale = qq))
}
if(s/500==floor(s/500)){
print(paste("s =",s, w, mu, sigma, logpost[s]))
}
}
## [1] "s = 500 0.65972914904991 4.28469003268875 0.405452035141612 -284.691001291928" "s = 500 0.34027085095009 1.99006962679139 0.250995183146609 -284.691001291928"
## [1] "s = 1000 0.595942847037826 4.27800551958405 0.377348480705164 -286.969076980901" "s = 1000 0.404057152962174 2.03894151043637 0.30982681142742 -286.969076980901"
## [1] "s = 1500 0.647284731246712 4.32792289508716 0.42383087208717 -283.984683957035" "s = 1500 0.352715268753288 2.06264939688377 0.293003644635685 -283.984683957035"
## [1] "s = 2000 0.681875974292433 4.31964666156271 0.411719777867089 -284.392728508546" "s = 2000 0.318124025707567 2.05510718985465 0.295371255576304 -284.392728508546"
## [1] "s = 2500 0.601103182730055 4.24400160743064 0.446830400212041 -285.459906307434" "s = 2500 0.398896817269945 2.03587612253421 0.241473527538476 -285.459906307434"
## [1] "s = 3000 0.661235314614079 4.28417774604423 0.399650525878003 -285.992706414487" "s = 3000 0.338764685385921 2.07827156671072 0.297205073536538 -285.992706414487"
## [1] "s = 3500 0.644527213106916 4.28812055333749 0.441452138427281 -283.120726285281" "s = 3500 0.355472786893084 2.01146104149005 0.258057420679898 -283.120726285281"
## [1] "s = 4000 0.687810647492561 4.2844578377926 0.458980223432612 -287.498421765689" "s = 4000 0.312189352507439 2.07489936081685 0.259843371397385 -287.498421765689"
## [1] "s = 4500 0.675398367857913 4.23091315341658 0.410829138479542 -285.115602530381" "s = 4500 0.324601632142087 2.00548178399042 0.263106146522126 -285.115602530381"
## [1] "s = 5000 0.623677983258687 4.2986070783669 0.460990184040201 -287.679382635602" "s = 5000 0.376322016741313 2.04273694828362 0.296539486159871 -287.679382635602"
## [1] "s = 5500 0.650641729955276 4.30297831389609 0.410588171771523 -284.091323207027" "s = 5500 0.349358270044724 2.05603614578286 0.314091813936972 -284.091323207027"
## [1] "s = 6000 0.691002769510701 4.30695193759895 0.43643062048752 -285.25059244604" "s = 6000 0.308997230489299 2.00724366068866 0.258042136715897 -285.25059244604"
## [1] "s = 6500 0.639653531193993 4.30228067466188 0.466022578492273 -289.071680298941" "s = 6500 0.360346468806007 2.0370347870673 0.223458414339796 -289.071680298941"
## [1] "s = 7000 0.657266286407271 4.25015359390783 0.408697090783594 -285.609121794953" "s = 7000 0.342733713592729 1.98352245342408 0.281123375854885 -285.609121794953"
## [1] "s = 7500 0.582150245729404 4.27752884035096 0.405335835898845 -288.834567817749" "s = 7500 0.417849754270596 2.03280535197236 0.302512079952417 -288.834567817749"
## [1] "s = 8000 0.597745067402946 4.30797847743645 0.426676105386858 -286.886008869559" "s = 8000 0.402254932597054 1.97410330661178 0.266695259035665 -286.886008869559"
## [1] "s = 8500 0.67048977917052 4.34923962430531 0.440934482727733 -287.976277284725" "s = 8500 0.32951022082948 1.99450886009298 0.253385375719033 -287.976277284725"
## [1] "s = 9000 0.676526152024463 4.32975227643208 0.424164941801625 -286.420409681512" "s = 9000 0.323473847975537 2.03912276672083 0.304347017719118 -286.420409681512"
## [1] "s = 9500 0.675252958080913 4.32408885533291 0.42428978538769 -284.557683174045" "s = 9500 0.324747041919087 2.01338154879308 0.269521567887484 -284.557683174045"
## [1] "s = 10000 0.657398814367743 4.2933157987876 0.458707092652969 -283.811172985985" "s = 10000 0.342601185632257 2.03282618126385 0.248515730931108 -283.811172985985"
## [1] "s = 10500 0.662058917060533 4.31483111471995 0.427552521301718 -284.108178852953" "s = 10500 0.337941082939467 2.00950256527503 0.279659447389632 -284.108178852953"
## [1] "s = 11000 0.669552336247041 4.36928639963658 0.427369327130037 -288.518332324442" "s = 11000 0.330447663752959 1.97992561769347 0.296908732685882 -288.518332324442"
## [1] "s = 11500 0.627693164995977 4.20937955935309 0.437031029444333 -287.287695169602" "s = 11500 0.372306835004023 2.0025029877848 0.248798427816352 -287.287695169602"
## [1] "s = 12000 0.706166055902923 4.32329970006521 0.423366944818915 -291.385728914099" "s = 12000 0.293833944097077 2.09369936916931 0.272464988230699 -291.385728914099"
Provide code to generate the density estimate on a grid consisting of 150 points in the interval [0,7].
Eruptions = seq(0,7, length=150)
## Compute the samples of the density over a dense grid
density.mcmc = array(0, dim=c(rrr-burn,length(Eruptions)))
for(s in 1:(rrr-burn)){
for(k in 1:KK){
density.mcmc[s,] = density.mcmc[s,] + w.out[s+burn,k]*dnorm(Eruptions,mu.out[s+burn,k],sigma.out[s+burn,k])
}
}
density.mcmc.m = apply(density.mcmc , 2, mean)
## Plot Bayesian estimate with pointwise credible bands along with kernel density estimate and frequentist point estimate
colscale = c("black", "blue", "red")
#yy = density(x)
yy = density(Eruptions)
density.mcmc.lq = apply(density.mcmc, 2, quantile, 0.025)
density.mcmc.uq = apply(density.mcmc, 2, quantile, 0.975)
plot(Eruptions, density.mcmc.m, type="n",ylim=c(0,max(density.mcmc.uq)),xlab="Duration", ylab="Density")
polygon(c(Eruptions,rev(Eruptions)), c(density.mcmc.lq, rev(density.mcmc.uq)), col="grey", border="grey")
lines(Eruptions, density.mcmc.m, col=colscale[3], lwd=2)
#lines(xx, density.EM, col=colscale[2], lty=2, lwd=2)
#lines(yy, col=colscale[3], lty=3, lwd=2)
points(x, rep(0,n))
Provide the a graph of the density estimate on the interval [0,7].
Is the sampler for \(c\) correct?
The full conditional for \(c_i\) is given by,
\[ \Pr(c_i = 1|...) = \frac{\omega\frac{1}{\sqrt[]{2\pi\sigma_1}}exp\left\{-\frac{1}{2\sigma^2_1}(x_i-\mu_1)^2\right\}}{\omega\frac{1}{\sqrt[]{2\pi\sigma_1}}exp\left\{-\frac{1}{2\sigma^2_1}(x_i-\mu_1)^2\right\}+(1-\omega)\frac{1}{\sqrt[]{2\pi\sigma_2}}exp\left\{-\frac{1}{2\sigma^2_2}(x_i-\mu_2)^2\right\}} \]
(This assumes that the mixture weights are parameterized as
for(i in 1:n){
v = rep(0,2)
v[1] = log(w) + dnorm(x[i], mu[1], sigma[1], log=TRUE) #Compute the log of the weights
v[2] = log(1-w) + dnorm(x[i], mu[2], sigma[2], log=TRUE) #Compute the log of the weights
v = exp(v - max(v))/sum(exp(v - max(v)))
cc[i] = sample(1:2, 1, replace=TRUE, prob=v)
}
Is the sampler for \(\omega\) correct?
The full conditional is simply
\[ \omega \mid \cdots \mbox{Beta} \left(\alpha_1 + \sum_{i=1}^{n} \mathbf{1}(c_i=1), \alpha_2 + \sum_{i=1}^{n} \mathbf{1}(c_i=2)\right) \]
The following one line of code implements it:
w = rbeta(1, aa[1] + sum(cc==1), aa[2] + sum(cc==2))
Are the sampler for \(\mu\) correct?
\[ \mu_k \mid \cdots \mbox{Normal}\left( \left[ \frac{n_k}{\sigma_k^2} + \frac{1}{\tau^2} \right]^{-1}\left[ \frac{\sum_{i:c_i=k} x_i}{\sigma_k^2} + \frac{\eta}{\tau^2} \right], \left[ \frac{n_k}{\sigma_k^2} + \frac{1}{\tau^2} \right]^{-1}\right) \]
The following code implements the sampler:
# Sample the means
for(k in 1:2){
nk = sum(cc==k)
xsumk = sum(x[cc==k])
tau2.hat = 1/(nk/sigma[k]^2 + 1/tau^2)
mu.hat = tau2.hat*(xsumk/sigma[k]^2 + eta/tau^2)
mu[k] = rnorm(1, mu.hat, sqrt(tau2.hat))
}
Remember, however, that there are various ways in which the implementation could proceed.
Are the samplers for \(\sigma_k\) correct?
The full conditional takes the form:
\[ 1/\sigma^2_k \mid \cdots \mbox{Gamma}\left(d + \frac{n_k}{2} , q + \frac{1}{2} \sum_{i:c_k=k} (x_i - \mu_k)^2 \right) \]
The following code implements the sampler:
# Sample the variances
for(k in 1:2){
nk = sum(cc==k)
dd.star = dd + nk/2
qq.star = qq + sum((x[cc==k] - mu[k])^2)/2
sigma[k] = sqrt(1/rgamma(1, dd.star, qq.star))
}
Remember, however, that there are various ways in which the implementation could proceed.
Does the code to generate the density estimate appear correct?
This is an example that, as before, assumes that the weights are parameterized as \(\omega\) and \(1- \omega\).
## Plot Bayesian estimate with pointwise credible bands
density.mcmc.lq = apply(density.mcmc, 2, quantile, 0.025)
density.mcmc.uq = apply(density.mcmc, 2, quantile, 0.975)
plot(xx, density.mcmc.m, type="n",ylim=c(0,max(density.mcmc.uq)),xlab="Eruptions", ylab="Density")
polygon(c(xx,rev(xx)), c(density.mcmc.lq, rev(density.mcmc.uq)), col="grey", border="grey")
lines(xx, density.mcmc.m, col="black", lwd=2)
points(x, rep(0,n))
Note that this will generate not only a point estimate, but also pointwise credible intervals.
Does the graph appear to be correct?
The density estimate should look similar to the Figure below. Note that the variance of both components is clearly different.
The full code for the MCMC algorithm follows:
x = faithful$eruptions
n = length(x)
plot(density(x))
points(x,rep(0,n))
## Priors set up using an "empirical Bayes" approach
aa = rep(1,2)
eta = mean(x)
tau = sqrt(var(x))
dd = 2
qq = var(x)/2
## Initialize the parameters
w = 1/2
mu = rnorm(2, mean(x), sd(x))
sigma = rep(sd(x)/2,2)
cc = sample(1:2, n, replace=T, prob=c(1/2,1/2))
## Number of iterations of the sampler
rrr = 12000
burn = 2000
## Storing the samples
cc.out = array(0, dim=c(rrr, n))
w.out = rep(0, rrr)
mu.out = array(0, dim=c(rrr, 2))
sigma.out = array(0, dim=c(rrr, 2))
logpost = rep(0, rrr)
for(s in 1:rrr){
# Sample the indicators
for(i in 1:n){
v = rep(0,2)
v[1] = log(w) + dnorm(x[i], mu[1], sigma[1], log=TRUE) #Compute the log of the weights
v[2] = log(1-w) + dnorm(x[i], mu[2], sigma[2], log=TRUE) #Compute the log of the weights
v = exp(v - max(v))/sum(exp(v - max(v)))
cc[i] = sample(1:2, 1, replace=TRUE, prob=v)
}
# Sample the weights
w = rbeta(1, aa[1] + sum(cc==1), aa[2] + sum(cc==2))
# Sample the parameters of the components
for(k in 1:2){
# Sample the means
nk = sum(cc==k)
xsumk = sum(x[cc==k])
tau2.hat = 1/(nk/sigma[k]^2 + 1/tau^2)
mu.hat = tau2.hat*(xsumk/sigma[k]^2 + eta/tau^2)
mu[k] = rnorm(1, mu.hat, sqrt(tau2.hat))
# Sample the variances
dd.star = dd + nk/2
qq.star = qq + sum((x[cc==k] - mu[k])^2)/2
sigma[k] = sqrt(1/rgamma(1, dd.star, qq.star))
}
# Store samples
cc.out[s,] = cc
w.out[s] = w
mu.out[s,] = mu
sigma.out[s,] = sigma
logpost[s] = 0
for(i in 1:n){
if(cc[i] == 1){
logpost[s] = logpost[s] + log(w) + dnorm(x[i], mu[1], sigma[1], log=TRUE)
}else{
logpost[s] = logpost[s] + log(1-w) + dnorm(x[i], mu[2], sigma[2], log=TRUE)
}
}
logpost[s] = logpost[s] + dbeta(w, aa[1], aa[2], log=TRUE)
for(k in 1:2){
logpost[s] = logpost[s] + dnorm(mu[k], eta, tau, log = T) + dgamma(1/sigma[k]^2, dd, qq)/sigma[k]^4
}
if(s/500==floor(s/500)){
print(paste("s =",s))
}
}
## [1] "s = 500"
## [1] "s = 1000"
## [1] "s = 1500"
## [1] "s = 2000"
## [1] "s = 2500"
## [1] "s = 3000"
## [1] "s = 3500"
## [1] "s = 4000"
## [1] "s = 4500"
## [1] "s = 5000"
## [1] "s = 5500"
## [1] "s = 6000"
## [1] "s = 6500"
## [1] "s = 7000"
## [1] "s = 7500"
## [1] "s = 8000"
## [1] "s = 8500"
## [1] "s = 9000"
## [1] "s = 9500"
## [1] "s = 10000"
## [1] "s = 10500"
## [1] "s = 11000"
## [1] "s = 11500"
## [1] "s = 12000"
## Compute the samples of the density over a dense grid
xx = seq(0,7,length=150)
density.mcmc = array(0, dim=c(rrr-burn,length(xx)))
for(s in 1:(rrr-burn)){
density.mcmc[s,] = density.mcmc[s,] + w.out[s+burn]*dnorm(xx,mu.out[s+burn,1],sigma.out[s+burn,1]) +
(1-w.out[s+burn])*dnorm(xx,mu.out[s+burn,2],sigma.out[s+burn,2])
}
density.mcmc.m = apply(density.mcmc , 2, mean)
## Plot Bayesian estimate with pointwise credible bands
density.mcmc.lq = apply(density.mcmc, 2, quantile, 0.025)
density.mcmc.uq = apply(density.mcmc, 2, quantile, 0.975)
plot(xx, density.mcmc.m, type="n",ylim=c(0,max(density.mcmc.uq)),xlab="Eruptions", ylab="Density")
polygon(c(xx,rev(xx)), c(density.mcmc.lq, rev(density.mcmc.uq)), col="grey", border="grey")
lines(xx, density.mcmc.m, col="black", lwd=2)
points(x, rep(0,n))
Is the sampler for \(c\) correct?
The full conditional for \(c_i\) is given by,
\[ \Pr(c_i = 1|...) = \frac{\omega\frac{1}{\sqrt[]{2\pi\sigma_1}}exp\left\{-\frac{1}{2\sigma^2_1}(x_i-\mu_1)^2\right\}}{\omega\frac{1}{\sqrt[]{2\pi\sigma_1}}exp\left\{-\frac{1}{2\sigma^2_1}(x_i-\mu_1)^2\right\}+(1-\omega)\frac{1}{\sqrt[]{2\pi\sigma_2}}exp\left\{-\frac{1}{2\sigma^2_2}(x_i-\mu_2)^2\right\}} \]
(This assumes that the mixture weights are parameterized as
for(i in 1:n){
v = rep(0,2)
v[1] = log(w) + dnorm(x[i], mu[1], sigma[1], log=TRUE) #Compute the log of the weights
v[2] = log(1-w) + dnorm(x[i], mu[2], sigma[2], log=TRUE) #Compute the log of the weights
v = exp(v - max(v))/sum(exp(v - max(v)))
cc[i] = sample(1:2, 1, replace=TRUE, prob=v)
}
Is the sampler for \(\omega\) correct?
The full conditional is simply
\[ \omega \mid \cdots \mbox{Beta} \left(\alpha_1 + \sum_{i=1}^{n} \mathbf{1}(c_i=1), \alpha_2 + \sum_{i=1}^{n} \mathbf{1}(c_i=2)\right) \]
The following one line of code implements it:
w = rbeta(1, aa[1] + sum(cc==1), aa[2] + sum(cc==2))
Are the sampler for \(\mu\) correct?
\[ \mu_k \mid \cdots \mbox{Normal}\left( \left[ \frac{n_k}{\sigma_k^2} + \frac{1}{\tau^2} \right]^{-1}\left[ \frac{\sum_{i:c_i=k} x_i}{\sigma_k^2} + \frac{\eta}{\tau^2} \right], \left[ \frac{n_k}{\sigma_k^2} + \frac{1}{\tau^2} \right]^{-1}\right) \]
The following code implements the sampler:
# Sample the means
for(k in 1:2){
nk = sum(cc==k)
xsumk = sum(x[cc==k])
tau2.hat = 1/(nk/sigma[k]^2 + 1/tau^2)
mu.hat = tau2.hat*(xsumk/sigma[k]^2 + eta/tau^2)
mu[k] = rnorm(1, mu.hat, sqrt(tau2.hat))
}
Remember, however, that there are various ways in which the implementation could proceed.
Are the samplers for \(\sigma_k\) correct?
The full conditional takes the form:
\[ 1/\sigma^2_k \mid \cdots \mbox{Gamma}\left(d + \frac{n_k}{2} , q + \frac{1}{2} \sum_{i:c_k=k} (x_i - \mu_k)^2 \right) \]
The following code implements the sampler:
# Sample the variances
for(k in 1:2){
nk = sum(cc==k)
dd.star = dd + nk/2
qq.star = qq + sum((x[cc==k] - mu[k])^2)/2
sigma[k] = sqrt(1/rgamma(1, dd.star, qq.star))
}
Remember, however, that there are various ways in which the implementation could proceed.
Does the code to generate the density estimate appear correct?
This is an example that, as before, assumes that the weights are parameterized as \(\omega\) and \(1- \omega\).
## Plot Bayesian estimate with pointwise credible bands
density.mcmc.lq = apply(density.mcmc, 2, quantile, 0.025)
density.mcmc.uq = apply(density.mcmc, 2, quantile, 0.975)
plot(xx, density.mcmc.m, type="n",ylim=c(0,max(density.mcmc.uq)),xlab="Eruptions", ylab="Density")
polygon(c(xx,rev(xx)), c(density.mcmc.lq, rev(density.mcmc.uq)), col="grey", border="grey")
lines(xx, density.mcmc.m, col="black", lwd=2)
points(x, rep(0,n))
Note that this will generate not only a point estimate, but also pointwise credible intervals.
Does the graph appear to be correct?
The density estimate should look similar to the Figure below. Note that the variance of both components is clearly different.
The full code for the MCMC algorithm follows:
x = faithful$eruptions
n = length(x)
plot(density(x))
points(x,rep(0,n))
## Priors set up using an "empirical Bayes" approach
aa = rep(1,2)
eta = mean(x)
tau = sqrt(var(x))
dd = 2
qq = var(x)/2
## Initialize the parameters
w = 1/2
mu = rnorm(2, mean(x), sd(x))
sigma = rep(sd(x)/2,2)
cc = sample(1:2, n, replace=T, prob=c(1/2,1/2))
## Number of iterations of the sampler
rrr = 12000
burn = 2000
## Storing the samples
cc.out = array(0, dim=c(rrr, n))
w.out = rep(0, rrr)
mu.out = array(0, dim=c(rrr, 2))
sigma.out = array(0, dim=c(rrr, 2))
logpost = rep(0, rrr)
for(s in 1:rrr){
# Sample the indicators
for(i in 1:n){
v = rep(0,2)
v[1] = log(w) + dnorm(x[i], mu[1], sigma[1], log=TRUE) #Compute the log of the weights
v[2] = log(1-w) + dnorm(x[i], mu[2], sigma[2], log=TRUE) #Compute the log of the weights
v = exp(v - max(v))/sum(exp(v - max(v)))
cc[i] = sample(1:2, 1, replace=TRUE, prob=v)
}
# Sample the weights
w = rbeta(1, aa[1] + sum(cc==1), aa[2] + sum(cc==2))
# Sample the parameters of the components
for(k in 1:2){
# Sample the means
nk = sum(cc==k)
xsumk = sum(x[cc==k])
tau2.hat = 1/(nk/sigma[k]^2 + 1/tau^2)
mu.hat = tau2.hat*(xsumk/sigma[k]^2 + eta/tau^2)
mu[k] = rnorm(1, mu.hat, sqrt(tau2.hat))
# Sample the variances
dd.star = dd + nk/2
qq.star = qq + sum((x[cc==k] - mu[k])^2)/2
sigma[k] = sqrt(1/rgamma(1, dd.star, qq.star))
}
# Store samples
cc.out[s,] = cc
w.out[s] = w
mu.out[s,] = mu
sigma.out[s,] = sigma
logpost[s] = 0
for(i in 1:n){
if(cc[i] == 1){
logpost[s] = logpost[s] + log(w) + dnorm(x[i], mu[1], sigma[1], log=TRUE)
}else{
logpost[s] = logpost[s] + log(1-w) + dnorm(x[i], mu[2], sigma[2], log=TRUE)
}
}
logpost[s] = logpost[s] + dbeta(w, aa[1], aa[2], log=TRUE)
for(k in 1:2){
logpost[s] = logpost[s] + dnorm(mu[k], eta, tau, log = T) + dgamma(1/sigma[k]^2, dd, qq)/sigma[k]^4
}
if(s/500==floor(s/500)){
print(paste("s =",s))
}
}
## [1] "s = 500"
## [1] "s = 1000"
## [1] "s = 1500"
## [1] "s = 2000"
## [1] "s = 2500"
## [1] "s = 3000"
## [1] "s = 3500"
## [1] "s = 4000"
## [1] "s = 4500"
## [1] "s = 5000"
## [1] "s = 5500"
## [1] "s = 6000"
## [1] "s = 6500"
## [1] "s = 7000"
## [1] "s = 7500"
## [1] "s = 8000"
## [1] "s = 8500"
## [1] "s = 9000"
## [1] "s = 9500"
## [1] "s = 10000"
## [1] "s = 10500"
## [1] "s = 11000"
## [1] "s = 11500"
## [1] "s = 12000"
## Compute the samples of the density over a dense grid
xx = seq(0,7,length=150)
density.mcmc = array(0, dim=c(rrr-burn,length(xx)))
for(s in 1:(rrr-burn)){
density.mcmc[s,] = density.mcmc[s,] + w.out[s+burn]*dnorm(xx,mu.out[s+burn,1],sigma.out[s+burn,1]) +
(1-w.out[s+burn])*dnorm(xx,mu.out[s+burn,2],sigma.out[s+burn,2])
}
density.mcmc.m = apply(density.mcmc , 2, mean)
## Plot Bayesian estimate with pointwise credible bands
density.mcmc.lq = apply(density.mcmc, 2, quantile, 0.025)
density.mcmc.uq = apply(density.mcmc, 2, quantile, 0.975)
plot(xx, density.mcmc.m, type="n",ylim=c(0,max(density.mcmc.uq)),xlab="Eruptions", ylab="Density")
polygon(c(xx,rev(xx)), c(density.mcmc.lq, rev(density.mcmc.uq)), col="grey", border="grey")
lines(xx, density.mcmc.m, col="black", lwd=2)
points(x, rep(0,n))
It’s useful to record some information about how your file was created.
suppressMessages(require('dplyr', quietly = TRUE))
suppressMessages(require('magrittr', quietly = TRUE))
suppressMessages(require('formattable', quietly = TRUE))
suppressMessages(require('knitr', quietly = TRUE))
suppressMessages(require('kableExtra', quietly = TRUE))
sys1 <- devtools::session_info()$platform %>%
unlist %>% data.frame(Category = names(.), session_info = .)
rownames(sys1) <- NULL
sys2 <- data.frame(Sys.info()) %>%
dplyr::mutate(Category = rownames(.)) %>% .[2:1]
names(sys2)[2] <- c('Sys.info')
rownames(sys2) <- NULL
if (nrow(sys1) == 9 & nrow(sys2) == 8) {
sys2 %<>% rbind(., data.frame(
Category = 'Current time',
Sys.info = paste(as.character(lubridate::now('Asia/Tokyo')), 'JST🗾')))
} else {
sys1 %<>% rbind(., data.frame(
Category = 'Current time',
session_info = paste(as.character(lubridate::now('Asia/Tokyo')), 'JST🗾')))
}
sys <- cbind(sys1, sys2) %>%
kbl(caption = 'Additional session information:') %>%
kable_styling(bootstrap_options = c('striped', 'hover', 'condensed', 'responsive')) %>%
row_spec(0, background = 'DimGrey', color = 'yellow') %>%
column_spec(1, background = 'CornflowerBlue', color = 'red') %>%
column_spec(2, background = 'grey', color = 'black') %>%
column_spec(3, background = 'CornflowerBlue', color = 'blue') %>%
column_spec(4, background = 'grey', color = 'white') %>%
row_spec(9, bold = T, color = 'yellow', background = '#D7261E')
rm(sys1, sys2)
sys
| Category | session_info | Category | Sys.info |
|---|---|---|---|
| version | R version 4.1.0 (2021-05-18) | sysname | Linux |
| os | Ubuntu 20.04.2 LTS | release | 5.8.0-54-generic |
| system | x86_64, linux-gnu | version | #61~20.04.1-Ubuntu SMP Thu May 13 00:05:49 UTC 2021 |
| ui | X11 | nodename | Scibrokes-Trading |
| language | en | machine | x86_64 |
| collate | en_US.UTF-8 | login | englianhu |
| ctype | en_US.UTF-8 | user | englianhu |
| tz | Asia/Tokyo | effective_user | englianhu |
| date | 2021-05-22 | Current time | 2021-05-22 04:27:25 JST🗾 |